% Copyright (C) 2009,2010,2011,2012  Marco Restelli
%
% This file is part of:
%   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
%
% LDGH is free software: you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% LDGH is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
% License for more details.
%
% You should have received a copy of the GNU General Public License
% along with LDGH. If not, see <http://www.gnu.org/licenses/>.
%
% author: Marco Restelli                   <marco.restelli@gmail.com>


function [p,e,t] = grid_gmsh2oct(d,gridfile,varargin)
% [p,e,t] = grid_gmsh2oct(d,gridfile,renumber)
%
% This function reads a .msh gmsh grid from file "gridfile" and
% reformats it in the usual [p,e,t] "Matlab-like" format.
%
% For 2D meshes, the third coordinate is ignored.
%
% renumber is optional and, if present, it defines a renumbering for
% the boundary labels in e. The length of renumber must be (at least)
% equal to the maximum boundary label used in gridfile. Consider that:
% a) setting renumber(i)<0 the sides with boundary label i are
%   discarded (this is useful, for instance, to drop internal
%   boundaries)
% b) if it's known that an intermediate label i is not used, the
%   corresponding element renumber(i) should be NaN, in order to
%   detect any error.
%
% LIMITATIONS: 
% - only one label for the domains is recognized
% - the curvilinear abscissa in e is not set

 fid = fopen(gridfile,'r');

 % Vertices
 fseek(fid,0,SEEK_SET);
 line = '';
 pos = [];
 while(isempty(pos))
   line = fgets(fid);
   pos = strfind(line,'$Nodes');
 end
 nv = fscanf(fid,'%d',1);
 % vertex coordinates
 xyz = fscanf(fid,'%d %f %f %f',[4,nv]);
 % It is possible that additional vertexes are included in the mesh
 % file which do not belong to any element. We thus check which are
 % the vertexes really belonging to the mesh while defining the
 % elements, and define p afterwards.
 vo_tmp = zeros(1,nv);
 vn_tmp = vo_tmp;
 nv = 0; % reset the counter

 % Elements and sides
 fseek(fid,0,SEEK_SET);
 line = '';
 pos = [];
 while(isempty(pos))
   line = fgets(fid);
   pos = strfind(line,'$Elements');
 end
 nes = fscanf(fid,'%d',1);
 line = fgets(fid); % skip an ampty line
 ns = 0;
 ne = 0;
 elpoint = ftell(fid);
 switch d % what are the elements and what are the sides?
  case 2
   telm = 2; % triangle
   tsde = 1; % line segment
  case 3
   telm = 4; % tetrahedron
   tsde = 2; % triangle
 endswitch
 for i=1:nes
   line = fgets(fid);
   ijk = str2num(line);
   switch ijk(2)
    case tsde
     ns = ns+1;
    case telm
     ne = ne+1;
    otherwise
     % nothing to do
   endswitch
 end
 t = zeros(d+2,ne);
 e = zeros(d+3,ns);
 ns = 0;
 ne = 0;
 fseek(fid,elpoint,SEEK_SET);
 for i=1:nes
   line = fgets(fid);
   ijk = str2num(line);
   switch ijk(2)
    case tsde
     ns = ns+1;
     % side vertexes
     e(1:d,ns) = ijk( 4+ijk(3) : 4+ijk(3)+(d-1) );
     e(d+3,ns) = ijk(    5     );
    case telm
     ne = ne+1;
     t(1:d+1,ne) = ijk( 4+ijk(3) : 4+ijk(3)+d );
     % renumber the vertexes
     for j=0:d
       jj = 4+ijk(3)+j;
       if(vo_tmp(ijk(jj))==0)
         nv = nv+1;
	 vo_tmp(ijk(jj)) = nv;
	 vn_tmp(nv) = ijk(jj);
       endif
     end
   endswitch
 end

 fclose(fid);

 % Now we can renumber the vertexes.
 % discard the label (and possibly the z coordinate)
 p = xyz(1+1:1+d,vn_tmp(1:nv));
 % use the new node order
 t(1:d+1,:) = vo_tmp(t(1:d+1,:));
 e(1:d  ,:) = vo_tmp(e(1:d  ,:));

 if(nargin>2)
   newl = varargin{1};
   side_del = [];
   for i=1:ns
     old_l = e(d+3,i);
     new_l = newl(old_l);
     if(new_l<0)
       side_del = [side_del i];
     else
       if(isnan(newl(e(d+3,i))))
         warning(['NaN label for side lable ',num2str(old_l)])
       endif
       e(d+3,i) = newl(e(d+3,i));
     end
   end
   e(:,side_del) = [];
 end

return

